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Abstract 


New  multiscale  methods  have  been  developed  for  the  study  of  microcracking  and 
dislocations.  The  multiscale  method  for  microcracking  aggregates  the  effect  of  an 
arbitrary  number  of  cracks  at  the  micro  level  into  a  single  equivalent  crack  at  the 
macrolevel.  The  methodology  is  based  on  the  concept  of  a  perforated  unit  cell,  which 
excludes  all  subdomains  of  the  unit  cell  where  the  material  loses  stability  in  the  averaging 
operations.  We  have  shown  that  for  this  approach,  the  macro  bulk  material  does  not  lose 
ellipticity  when  failure  occurs.  This  is  an  important  development  since  it  avoids  the  ill- 
posedness  usually  associated  with  multiscale  modeling  of  failure.  Comparisons  with 
direct  numerical  simulations  show  excellent  accuracy.  We  have  also  developed  new 
models  for  dislocations  and  concurrent  methods  that  couple  atomistic  models  with  these 
new  continuum  dislocation  methods.  These  enable  dislocation  dynamics  in  arbitrary 
geometries  and  materials  to  be  modeled  in  terms  of  fundamental  physics. 


Statement  of  the  Problem  Studied 

Two  problems  were  studied  in  this  investigation: 

1 .  the  multiscale  modeling  of  failure 

2.  the  multiscale  modeling  of  dislocations 

The  multiscale  modeling  of  failure  remains  an  unresolved  problem  because  when  a 
subscale  model  predicts  that  the  stress-strain  behavior  at  the  coarser  scale  loses  ellipticity, 
the  coarse  scale  model  becomes  ill-posed.  Consequently,  the  results  of  the  coarse-scale 
model  become  very  mesh  dependent  and  unreliable.  This  occurs  for  example  when 
microcracking  takes  place  at  the  finer  scale.  The  objective  of  this  investigation  was  to 
develop  methods  for  circumventing  these  difficulties  by  separating  out  the  unstable 
material  behavior  at  the  finer  scale  into  an  equivalent  discontinuity  and  traction.  The 
failure  behavior  at  the  fine  scale  is  thus  isolated  as  either  a  crack,  shear  band  or 
dislocation  at  the  coarser  scale.  This  is  in  agreement  with  observed  physics  of  failure 
processes,  since  they  are  always  associated  with  the  aforementioned  phenomena. 

Dislocations  are  closely  related  phenomena  which  are  associated  with  plasticity. 
Dislocations  are  characterized  by  relative  atomic  displacements  at  the  crystalline  level 
which  can  be  characterized  by  discontinuities  at  the  macrolevel.  In  this  work,  new 
methods  have  been  developed  for  modeling  dislocations  by  the  extended  finite  element 
method.  Furthermore,  methods  have  been  developed  for  coupling  continuum  models  of 
dislocations  with  atomistic  models.  These  coupled  multiscale  models  will  allow 
investigators  to  bypass  the  heuristic,  phenomenological  rules  which  are  used  to  drive  the 
displacements  of  dislocation  cores  in  dislocation  dynamics  models. 


Summary  of  Important  Results 

It  has  been  found  that  ill-posedness  can  be  avoided  in  both  hierarchical  and 
semiconcurrent  couplings  of  multiscale  models  through  the  use  of  two  key  concepts.  The 
first  key  concept  is  the  definition  of  a  perforated  representative  volume  element  (or  unit 
cell)  which  excludes  all  material  points  that  are  unstable  (to  be  more  precise,  not  convex). 
The  second  key  concept  is  a  method  for  defining  an  equivalent  discontinuity  in  terms  of 
the  behavior  of  the  perforated  unit  cell.  The  method  is  called  the  multiscale  aggregating 
discontinuity  method  (MAD). 

Let  the  perforated  subdomain  of  the  unit  cell  Q”  be  denoted  by  Q”  .  The  macroscale  first 

Piola-Kirchhoff  stress  is  then  defined  in  terms  of  the  microscale  first  Piola-Kirchhoff 
stress  P”  by 
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Similarly  the  macroscale  deformation  gradient  is  defined  in  terms  of  the  microscale 
deformation  gradient  F”  by 
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Let  the  motion  (jf  (X™ ,  t)  on  the  sides  of  the  unit  cell  be  given  by 

=  onX"’edn; 


(1.3) 


We  have  shown  that  an  equivalent  discontinuity  at  the  macrolevel  U'^  with  normal 
can  then  be  extracted  by 
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This  macroscale  discontinuity  is  then  injected  into  the  macroscale  model.  The  forces 
required  to  enforce  this  discontinuity  are  the  traction  across  the  macrocrack. 

We  have  shown  that  the  macroscale  model  is  elliptic  as  long  if  the  perforated  unit  cell 
excludes  all  material  that  is  not  convex.  To  put  this  in  more  precise  terms,  recall  that  the 
convexity  of  the  micromaterial  requires  that 


VF” 


F”  :  C”  :  F”  >  0 


(1.5) 


where  C”  is  the  tangent  modulus  for  the  micromaterial.  Ellipticity  of  the  macromaterial 
requires  that 


g0h:C^:g0h>O  Vg  and  h 


(1.6) 


It  has  been  shown  in  Belytschko  et  al.  [1]  that  Equation  (1.5)  implies  (1.6)  for  the  method 
as  described.  Therefore  the  multiscale  aggregating  equivalent  discontinuity  method 
guarantees  well-posedness. 

The  evaluation  of  a  multiscale  method  for  failure  poses  some  particular  challenges. 
Comparison  with  experimental  results  is  usually  inconclusive,  since  the  material 
constants  for  failure  are  not  readily  available.  We  were  not  able  to  find  experiments  that 
provide  both  subscale  properties  and  a  macroscale  response.  Therefore,  we  have  chosen 
to  use  comparison  of  the  multiscale  results  with  direct  numerical  simulation  (DNS)  to 
evaluate  the  method.  We  show  one  example  here:  a  four-point  bending  problem  of  a 
composite  beam. 

The  problem  is  shown  in  Figure  1.  The  fibers  are  oriented  normal  to  the  plane  of  the 
beam.  The  beam  is  initially  notched  at  the  center  on  the  bottom  surface  as  shown.  Figure 
1  also  shows  the  unit  cells  for  the  multiscale  model.  Each  unit  cell  is  represented  by  a 
finite  element  in  the  coarse  scale  model. 
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Figure  1.  Initial  setup  for  failure  of  composite  beam  due  to  four-point  bending, 


Figure  2(a)  shows  the  DNS  model,  which  represents  the  detailed  structure  of  the 
composite.  The  crack  path  in  the  DNS  model  is  shown  in  Figure  2(b).  Note  that  the  crack 
path  for  the  DNS  is  quite  jagged,  since  it  bypasses  many  of  the  fibers.  The  coarse-grained 
solution  is  shown  in  Figure  2(c).  Here  the  crack  path  is  much  smoother,  because  the 
subscale  structure  in  the  unit  cells  is  not  apparent  to  the  coarse  scale  computation. 
However,  the  overall  crack  path  is  quite  similar. 


(C) 

Figure  2.  Four-point  bending  beam  probiem:  (a)  flnite  element  discretization  for  the  DNS,  (b)  the 
plots  of  effective  stress  contour  for  the  DNS,  and  (c)  the  plots  of  effective  stress  contour  for  the 

coarse-grained  MAD  method. 

The  adequacy  of  the  coarse  scale  computation  by  our  method  can  be  judged  by  the  force- 
deflection  curve  shown  in  Figure  3.  We  first  compare  the  MAD  result  with  the  hourglass 
mode  (which  is  a  recent  improvement  in  the  MAD  method  that  accounts  for  hourglassing 
in  the  unit  cell,  see  [1]).  In  that  case,  the  agreement  with  the  DNS  solution  is 
exceptionally  good:  although  there  are  some  deviations  in  the  time  history  that  are  as 
large  as  15%,  the  peak  load  is  reproduced  very  well.  The  older  method,  which  does  not 
account  for  the  hourglass  deformations,  agrees  better  in  the  early  part  of  the  response,  but 
the  force  overshoots  the  DNS  force  significantly,  which  indicates  that  without  the 
hourglass  deformation  modes  the  method  is  too  stiff. 

As  can  be  seen,  the  multiscale  method  reproduces  the  fine  scale  result  quite  accurately.  It 
does  so  with  a  speedup  of  approximately  30.  This  is  the  major  benefit  of  multiscale 
methods,  for  without  such  speedups,  larger  problems  are  beyond  the  reach  of  even  the 
largest  computers. 


Figure  3.  Comparison  of  load-deflection  curves  between  the  DNS  and  the  MAD  for  the  four-point 

bending  beam  problem. 


The  new  method  for  modeling  dislocations  is  based  on  introducing  discontinuities  into 
the  displacement  field  in  the  spirit  of  Volterra’s  concept  of  dislocations  and  the  extended 
finite  element  method  (XFEM)  framework.  The  displacement  field  surrounding  a 
network  of  dislocations  in  a  material  is  represented  in  a  finite  element  model  by  the 
combination  of  standard  piecewise  continuous  polynomial  basis  functions  and 
discontinuous  enrichment  functions.  The  enrichment  functions  are  nonzero  only  in 
elements  which  have  been  intersected  by  the  slip  plane,  and  since  the  relative  slip  across 
the  glide  plane  is  fixed  to  a  translation  vector  of  the  lattice,  the  enrichment  adds  no 
additional  degrees  of  freedom  to  the  discretized  finite  element  system  of  equations  but 
instead  can  be  treated  as  an  additional  set  of  nodal  forces.  The  advantage  of  this 
approach  is  that  it  can  be  used  with  complex  geometries,  anisotropic  and  non-uniform 
materials,  and  nonlinear  materials  [2]. 

A  key  feature  of  this  method  is  the  use  of  level  set  functions  to  describe  the  location  of 
slip  surfaces  and  dislocation  lines,  and  the  enrichment  functions  that  inject  the 
discontinuity  across  the  glide  plane  to  the  solution.  Each  dislocation  is  defined  by  two 
scalar  fields;  the  first, /(xj,  is  a  signed  distance  function  to  the  slip  surface,  and  the  second, 
g(x),  gives  the  signed  distance  to  the  dislocation  line  from  the  projection  of  every 
material  point  to  the  slip  surface.  Figure  4  shows  how,  for  a  curved  dislocation  line,  the 
level  set  functions  /  and  g  vary  and  make  up  a  local  coordinate  system  around  the 
dislocation  core. 


wt 


Figure  4.  Illustration  of  the  local  coordinate  system  deflned  about  the  core  of  a  dislocation  by  the 

level  set  functions  f(x)  and  g(x). 

The  second  major  component  of  this  new  method  for  dislocations  is  the  construction  of 
the  enrichment  functions.  The  simplest  possibility  is  a  constant  jump  across  the  glide 
surface  following  the  Volterra  construction.  In  addition  to  its  simplicity,  since  the 
function  is  constant  everywhere  except  for  the  discontinuity,  there  are  no  enrichment 
strains  to  be  integrated.  This  is  especially  helpful  when  the  dislocation  line  is  curved 
because  it  eliminates  the  enrichment  strain  caused  by  the  rotation  of  the  local  level  set 
coordinate  system.  As  shown  in  Oswald  et  al.  [3]  a  more  sophisticated  option  for 
enrichment  is  to  take  the  infinite  domain  solutions  to  dislocations  near  the  core  using  the 
local  coordinate  system  defined  by  the  level  set  functions.  The  analytical  solutions  can 
dramatically  improve  the  accuracy  of  the  computation,  as  shown  in  Figure  5  where  a 
27,000  element  solution  with  singular  enrichments  has  comparable  accuracy  to  a  two- 
million  element  solution  with  the  jump  enrichment.  Both  of  the  solutions  employed  a 
uniform  mesh,  so  with  an  adaptively  generated  mesh  the  number  of  degrees  of  freedom 
for  the  jump  enrichment  could  be  significantly  reduced. 
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Figure  5.  Comparison  of  Oxx  stress  values  on  a  line  parallel  to  the  x-axis  passing  above  a  circular 
dislocation  loop  lying  in  the  x-y  plane  of  an  isotropic  material  for  a  constant  jump  enrichment  and 
singular  enrichments  formed  from  the  solutions  of  infinite-length  straight  dislocations. 


As  an  illustration  of  the  application  of  this  approach,  consider  the  silicon-germanium 
quantum  dot  shown  in  Figure  6.  This  problem  features  several  of  the  difficulties  faced  by 
current  dislocation  modeling  methods.  The  quantum  dot  geometry  considered  here  is 
composed  of  36  facets  representing  9  unique  energetically  favorable  crystalline  surfaces. 


Both  the  quantum  dot  and  the  silicon  substrate  are  single  crystals,  and  therefore  exhibit 
considerable  anisotropic  effects.  Additionally,  as  reported  by  Malachias  et  al.  [4]  the  Ge 
concentration  varies  from  nearly  pure  germanium  at  the  top  of  the  quantum  dot  to  a 
silicon  rich  core  at  the  base.  As  a  result,  not  only  is  there  a  material  interface  to  the 
substrate  at  the  base,  but  the  elastic  properties  vary  spatially  throughout  the  quantum  dot. 
Previous  elastic  analyses  of  dislocations  in  these  structures  have  required  the  assumption 
of  uniform  isotropic  material  properties,  but  the  XFEM  method  can  treat  these  issues 
without  complication. 


Figure  6.  Finite  element  model  of  a  dislocation  loop  at  the  interface  of  a  SiGe  quantum  dot  deposited 
on  a  (001)  Si  surface  (left),  and  the  resulting  magnitude  of  the  stress  generated  by  the  lattice 
mismatch  between  SiGe  and  Si,  and  distortions  in  the  lattice  caused  by  the  dislocation  loop. 


Another  way  that  this  dislocation  method  can  be  extended  is  by  modeling  the  core 
structure  of  the  dislocation  by  an  atomistic  method  such  as  molecular  statics  or  molecular 
dynamics.  The  modifications  needed  to  achieve  this  are  a  suitable  physical  model  for  the 
interactions  between  atoms,  a  nonlinear  minimization  solver  for  molecular  statics  or  a 
time  integrator  for  molecular  dynamics,  and  a  coupling  method  to  the  continuum 
description  of  the  material.  To  demonstrate  this  extension,  a  dislocation  dipole  ejected  by 
a  vacancy  was  modeled  in  a  graphene  sheet,  where  the  carbon-carbon  interaction  energies 
are  approximated  by  the  Tersoff-Brenner  potential.  A  modified  version  of  the  bridging 
domain  method  (BDM)  enforces  compatibility  between  the  continuum  and  atomistic 
regions.  The  modification  to  the  BDM  allows  the  discontinuity  to  intersect  the  coupling 
domain  such  that  it  perfectly  matches  the  discontinuity  of  atomic  displacements  along  the 
slipped  edges  of  the  graphene.  Figure  7  shows  the  computed  displacement  fields  for  the 
dislocation  dipole  and  vacancy,  and  the  XFEM  enriched  elements  encompass  the  glide 
plane  in  the  continuum. 


This  new  method  for  dislocations  makes  possible  the  modeling  of  dislocations  in  very 
complex  geometries  and  arbitrary  material,  including  both  anisotropy  and  nonlinearities. 
Furthermore,  it  is  ideally  suited  for  combined  atomistic,  continuum  solutions,  so  it  can 
efficiently  solve  large  scale  problems  where  first  principles  models  are  used  for  the  core. 
This  will  avoid  the  need  for  phenomenological  rules  for  the  motion  of  the  dislocations. 
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Figure  7.  Modeling  of  a  dislocation  ejected  by  a  vacancy  in  graphene  sheet  by  the  extended  flnite 
element  method  conpled  to  a  molecular  mechanics  simnlation  in  the  defect  regions. 
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